/*
energyTable:  brief Definition of EnergyTable class.
    Copyright (C) 2010  Jorge E. A. Atala
                        Adrian Remonda
                        FuDePAN
    This file is part of the Backbones-Selector Proyect.

    Backbones-Selector is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Backbones-Selector is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Backbones-Selector.  If not, see <http://www.gnu.org/licenses/>.
*/


/*
layout del archivo de aminoacidos
cantidad de aminoacidos
aminoacido 1
aminoacido 2
...
aminoacido m
cantidad de rangos de distancias
extremo inferior del primer rango
extremo superior del primer rango
extremo superior del segundo rango
...
extremo superior de rango n
datos de energia de la primera matriz en formato triangular
datos de energia de la segunda matriz en formato triangular
...
datos de energia de la enesima matriz en formato triangular
*/

#ifndef ENERGYTABLE_H
#define ENERGYTABLE_H

#include <vector>
#include <string>
#include <algorithm>
#include <iostream>
#ifndef MILI_NAMESPACE
#   define MILI_NAMESPACE
#endif
#include <mili/mili.h>
#include <fstream>
//#define NO_BIO_MOLECULAR
#define NO_RNA_STRUCTURE
#include <biopp/biopp.h>
using namespace biopp;


namespace backbonesSelector
{


//------------------------------------------------------------------------------------
//Begin of EnergyTable Class
//------------------------------------------------------------------------------------
const unsigned int AMINOACID_TO_MATRIX_TABLE_SIZE = 256;

template < class Energy = float, class Distance = float >
class EnergyTable
{
private:
    typedef unsigned int MatrixIndex;
    Energy*** energyMatrix;
    std::vector<AminoacidStd> aminoacidsList;
    size_t aminoacidsCount;
    std::vector<Distance> rangesList;
    size_t rangesCount;
    inline MatrixIndex aminoAcidToMatrixPosition(const char a) const;
    void buildAminoacidToMatrixTable();
    struct sAminoacidToMatrixItem
    {
        bool present;
        MatrixIndex index;
        sAminoacidToMatrixItem(): present(false), index(0) {}
    };
    std::vector<sAminoacidToMatrixItem> aminoacidToMatrixTable;
public:
    EnergyTable(const char* fileData);
    inline Distance getMaxDistance() const;
    Energy getEnergy(const char aminoacid1, const char aminoacid2, const Distance distance) const;
    ~EnergyTable();
};

//Constructor of EnergyTable class
//parameter "ranges": collection of upper bounds for each range of distances
template<class Energy, class Distance>
EnergyTable<Energy, Distance>::EnergyTable(const char* fileData):
    aminoacidToMatrixTable(AMINOACID_TO_MATRIX_TABLE_SIZE)
{
    std::ifstream f;
    f.exceptions(std::ifstream::failbit | std::ifstream::badbit | std::ifstream::eofbit);

    try
    {
        f.open(fileData);

        std::cout << "Opening the tables file successfully" << std::endl;

        //read aminoacids
        f >> aminoacidsCount;

        aminoacidsList.resize(aminoacidsCount);
        for (unsigned int it = 0; it < aminoacidsCount; ++it)
        {
            char temp;
            f >> temp;
            aminoacidsList[it] = temp;
        }

        //the table max value is the highest aminoacid value
        buildAminoacidToMatrixTable();

        //read ranges distances
        f >> rangesCount;

        rangesList.resize(rangesCount + 1);
        for (unsigned int it = 0; it < rangesCount + 1; ++it)
            f >> rangesList[it];

        //read the energy matrix
        unsigned int i;
        energyMatrix = new Energy** [rangesCount];
        for (i = 0; i < rangesCount; ++i)
        {
            energyMatrix[i] = new Energy*[aminoacidsCount];
            for (unsigned int x = 0; x < aminoacidsCount ; ++x)
                energyMatrix[i][x] = new Energy[aminoacidsCount];
        }

        for (i = 0; i < rangesCount; ++i)
            for (unsigned int x = 0; x < aminoacidsCount; ++x)
                for (unsigned int y = x; y < aminoacidsCount; ++y)
                {
                    f >> energyMatrix[i][x][y];
                    energyMatrix[i][y][x] = energyMatrix[i][x][y]; //when i do a triangular matrix, this will be not necesary
                }
    }
    catch (const std::ifstream::failure& e)
    {
        std::cerr << "Exception opening/reading energy tables file" << std::endl;
        std::cerr << e.what() << std::endl;
    }
}

template<class Energy, class Distance>
EnergyTable<Energy, Distance>::~EnergyTable()
{
    for (unsigned int i = 0; i < rangesCount; ++i)
    {
        for (unsigned int x = 0; x < aminoacidsCount; ++x)
            delete [] energyMatrix[i][x];
        delete [] energyMatrix[i];
    }
    delete [] energyMatrix;
}

template<class Energy, class Distance>
void EnergyTable<Energy, Distance>::buildAminoacidToMatrixTable()
{
    MatrixIndex index = 0;

    for (MatrixIndex i = 0; i < aminoacidsCount; ++i)
    {
        index = aminoacidsList[i].as_char();
        aminoacidToMatrixTable[index].present = true;
        aminoacidToMatrixTable[index].index = i;
    }
}

//This function gets the aminoacid name and return the corresponding row/column position into the energyMatrix
template<class Energy, class Distance>
inline typename EnergyTable<Energy, Distance>::MatrixIndex EnergyTable<Energy, Distance>::aminoAcidToMatrixPosition(const char a)const
{
    if (!aminoacidToMatrixTable[MatrixIndex(a)].present)
        throw AminoacidNotFound();

    return aminoacidToMatrixTable[MatrixIndex(a)].index;
}

//Obtain the energy value from the energyMatrix, corresponding whit the parameters values
template<class Energy, class Distance>
Energy EnergyTable<Energy, Distance>::getEnergy(const char aminoacid1, const char aminoacid2, const Distance distance) const
{
    Energy result = 0;

    const unsigned int energyRow = this->aminoAcidToMatrixPosition(aminoacid1);
    const unsigned int energyColumn = this->aminoAcidToMatrixPosition(aminoacid2);

    if (mili::in_range(distance, rangesList.front(), rangesList.back()))
    {
        typename std::vector<Distance>::const_iterator boundPosition;
        //The depth of the cell in the matrix is obtained with the distance between the aminoacids,
        //we need the upper bound position in "distanceRanges"
        boundPosition = lower_bound(rangesList.begin(), rangesList.end(), distance);
        unsigned int energyDepth = boundPosition - rangesList.begin();

        //If we have an energyDepth equal to minimal distance of first range,
        //we need plus 1 for put this param into the first range, because the
        //first element inth array is the minimal distance for aminoacids interaction,
        //not like the other elements, they are the max distance of the previous
        //range too.
        if (energyDepth == 0)
            energyDepth = 1;
        result = energyMatrix[energyDepth - 1][energyRow][energyColumn];
    }
    return result;
}

//Returns the max aminoacids interaction distance
template<class Energy, class Distance>
Distance EnergyTable<Energy, Distance>::getMaxDistance() const
{
    return rangesList.back();
}
}

//------------------------------------------------------------------------------------
//End of EnergyTable Class
//------------------------------------------------------------------------------------
#endif
